Load packages and configure the grid resolution and source filename:
library(dplyr)
library(ggplot2)
library(dggridR)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(prettygrids)
world <- ne_countries(scale = "medium", returnclass = "sf")
resolution <- 7
filename <- "occurrence_minimal_20200212"
Use the dggridR package to create a discrete global grid, and fetch the grid cell surface area to be included in the map title:
dggs <- dgconstruct(projection = "ISEA", res = resolution, resround = "up")
area <- format(signif(dggetres(dggs)$area_km[resolution + 1], 2), big.mark = ",")
This notebook uses a A CSV file with occurrence records exported from the database database. Alternatively the robis R package can be used, but downloading the entire database can take a while. This downloads the zipped CSV file and aggregates the data into locations with number of records:
if (!file.exists("grouped.temp")) {
if (!file.exists("data.zip")) {
download.file(paste0("https://download.obis.org/export/", filename, ".zip"), "data.zip")
}
df_raw <- read.csv(unz("data.zip", paste0(filename, ".csv")), stringsAsFactors = FALSE)
df_grouped <- df_raw %>%
group_by(decimallongitude = round(decimallongitude, 2), decimallatitude = round(decimallatitude, 2), speciesid) %>%
summarize(records = n(), species = length(unique(speciesid)))
save(df_grouped, file = "grouped.temp")
file.remove("data.zip")
} else {
load("grouped.temp")
}
Next we assign cell IDs to each location and caluclate statistics per cell:
df_grouped$cell <- dgtransform(dggs, df_grouped$decimallatitude, df_grouped$decimallongitude)
stats <- df_grouped %>%
group_by(cell) %>%
summarise(records = sum(records), species = length(unique(speciesid)))
map_theme <- theme(
axis.ticks.length = unit(0, "pt"),
line = element_blank(), rect = element_blank(),
axis.text = element_blank(), axis.title = element_blank()
)
map_scale_records <- scale_fill_viridis_c(
option = "plasma",
trans = "log10",
name = "Records",
breaks = c(1, 10, 100, 1000, 10000),
labels = function(x) format(x, scientific = FALSE),
limits = c(1, 100000),
oob = scales::squish
)
map_scale_species <- scale_fill_viridis_c(
option = "plasma",
trans = "log10",
name = "Species",
labels = function(x) format(x, scientific = FALSE)
)
Depending on the projection the grids produced by dggridR can cause artefacts in your map. Here we are using the prettygrids::make_grid utility function (still under development) to fix some of these problems.
grid_pac <- make_grid(offset = -180, res = resolution) %>%
merge(stats, by.x = "cell", by.y = "cell")
crs <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
bbox <- st_bbox(st_as_sfc(st_bbox(c(xmin = 18, xmax = 248, ymin = -55, ymax = 55), crs = 4326)) %>% st_transform(crs))
ggplot() +
geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = records), size = 0, color = NA) +
geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs, xlim = c(bbox$xmin, bbox$xmax), ylim = c(bbox$ymin, bbox$ymax)) +
map_theme +
map_scale_records +
labs(
title = "Number of records",
subtitle = paste0("Number of records per grid cell of approximately ", area, " km2, Pacific centered Robinson projection. Ocean Biogeographic Information System (OBIS), 2020.")
)
ggplot() +
geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = species), size = 0, color = NA) +
geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs, xlim = c(bbox$xmin, bbox$xmax), ylim = c(bbox$ymin, bbox$ymax)) +
map_theme +
map_scale_species +
labs(
title = "Number of species",
subtitle = paste0("Number of species per grid cell of approximately ", area, " km2, Pacific centered Robinson projection. Ocean Biogeographic Information System (OBIS), 2020.")
)
north <- st_set_crs(st_as_sf(as(raster::extent(-180, 180, 0, 90), "SpatialPolygons")), 4326)
grid_ortho <- make_grid(offset = NA, res = resolution) %>%
merge(stats, by.x = "cell", by.y = "cell") %>%
filter(st_intersects(geometry, north, sparse = FALSE))
crs <- "+proj=laea +lat_0=52 +lon_0=0 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_defs"
bbox <- st_sfc(st_polygon(list(matrix(c(-140, 40, -60, 40, 60, 40, 140, 40, -140, 40), ncol = 2, byrow = TRUE))), crs = 4326)
bbox_trans <- bbox %>% st_transform(crs)
ggplot() +
geom_sf(data = grid_ortho %>% st_transform(crs), aes(fill = records), size = 0, color = NA) +
geom_sf(data = world %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs,
xlim = c(min(st_coordinates(bbox_trans)[,1]), max(st_coordinates(bbox_trans)[,1])),
ylim = c(min(st_coordinates(bbox_trans)[,2]), max(st_coordinates(bbox_trans)[,2]))) +
map_theme +
map_scale_records +
labs(
title = "Number of records",
subtitle = paste0("Number of records per grid cell of approximately ", area, " km2, Lambert azimuthal equal-area projection. Ocean Biogeographic Information System (OBIS), 2020.")
)
ggplot() +
geom_sf(data = grid_ortho %>% st_transform(crs), aes(fill = species), size = 0, color = NA) +
geom_sf(data = world %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
coord_sf(crs = crs,
xlim = c(min(st_coordinates(bbox_trans)[,1]), max(st_coordinates(bbox_trans)[,1])),
ylim = c(min(st_coordinates(bbox_trans)[,2]), max(st_coordinates(bbox_trans)[,2]))) +
map_theme +
map_scale_species +
labs(
title = "Number of species",
subtitle = paste0("Number of species per grid cell of approximately ", area, " km2, Lambert azimuthal equal-area projection. Ocean Biogeographic Information System (OBIS), 2020.")
)